package hr.irb.goTermList;

import hr.irb.geneOntology.GOTerm;
import hr.irb.geneOntology.GeneOntology;
import hr.irb.geneOntology.GoNamespace;
import hr.irb.geneOntology.GoTermProperties;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.LinkedHashMap;
import java.util.List;
import java.util.Locale;
import java.util.Map;

/**
 * Runs Revigo from the command line on a single file, or many files in a given 
 * directory.
 * 
 * @author Fran Supek
 */
public class BatchRevigoRunner {
    
  /** Cutoff for the pairwise semantic similarity. 
   * Smaller number = shorter list. */
  public static double cutoff = 0.7; 
  
  /** JIANG, LIN, RESNIK or SIMREL. */
  public static String measure = "SIMREL";
  
  /**
   * Can be "p_value", "higher", "lower", "absolute" or "abs_log".
   */
  public static String whatIsBetter = "p_value";
  
  public static final String goaPackageName = "goa_2014_09_30";
  
  
  private static int goIdColumnIdx = 1;
  private static int valueColumnIdx = 2;
  
  /** The 9606 is human, the 0 is whole Uniprot. See revigo.jsp for others. */
  public static int taxId = 0;
  
  /** If true, prints to stdout and not to a file. */
  public static boolean stdout = false;
  
  
  /**
   * Processes a file or a directory using REVIGO and writes 
   * the output files to the same directory. 
   */
  public static void main(String[] args) throws IOException, ClassNotFoundException, Exception {
  
    if ( args.length==0 || args[0].startsWith("-") ) {
      // prints usage instructions
      System.out.println(
"java -jar revigo.jar fileOrDirToProcess [--goIdColumnIdx=n] \n" +
"[--valueColumnIdx=n] [--cutoff=n] [--taxId=n] " + 
"[--whatIsBetter=p_value|higher|lower|absolute|abs_log] " +
"[--measure=SIMREL|RESNIK|LIN|JIANG] [--stdout] \n" +
"\n" +
"fileOrDirToProcess Tab-separated text table(s). Any rows starting with \"#\" \n" +
"or \"%\" are ignored. First row is not ignored, unless starting with one of the\n" +
"given characters. For each input file, creates one output file with the\n" +
"appended suffix '_revigoOut.txt'\n" +
"\n" +
"[--goIdColumnIdx=n] Index of the column with the GO IDs. Default = first\n" +
"column (1)\n" +
"\n" +
"[--valueColumnIdx=n] Index of the column with the p-values, enrichments or\n" +
"similar. Default = second column (2). Note: all other columns in the table are\n" + 
"ignored. Pass 0 to work without p-values (not recommended).\n" +
"\n" +
"[--cutoff=n] Cutoff for the pairwise semantic similarity. Smaller values =\n" + 
"shorter list. Settings below 0.5 tend to be too restrictive and are not\n" + 
"recommended. Default = 0.7.\n" +
"\n" +
"[--taxId=n] The database of GO terms' frequencies to be used in calculating\n" +
"semantic similarity. Default = 0 (uses whole Uniprot). This is the NCBI\n" +
"taxonomy ID: E. coli is 83333, human is 9606, etc. For a list of the supported\n" +
"species, refer to directory /hr/irb/goTermList/goTermSizesData/goa_2013_04_30\n" +
"in the JAR archive. If your species is absent, it's generally safe to use 0.\n" +
"\n" +
"[--whatIsBetter=p_value|higher|lower|absolute|abs_log] How the data in the\n" +
"\"valueColumn\" is treated. Default is \"p_value\". For enrichments or similar,\n" +
"use \"abs_log\". For other quantities, choose whether lower, higher, or higher\n" +
"absolute values are more desirable.\n" +
"\n" +
"[--measure=SIMREL|RESNIK|LIN|JIANG] The semantic similarity measure to use.\n" +
"Default = SIMREL.\n" +
"\n" +
"[--stdout] Prints to standard output instead of to a file (default behavior).\n"        
      );
      return;
    }
    
    String inputDir = args[0];
    for (int i=1; i<args.length; i++){
      String arg = args[i].toLowerCase();
      String argName = arg.split("=")[0];
      String argVal = ( arg.split("=").length == 1 ? null : arg.split("=")[1]);
      if ( argName.equals("--taxid") ) {
        
      } else if ( argName.equals("--goidcolumnidx") ) {
        goIdColumnIdx = Integer.parseInt(argVal);
      } else if ( argName.equals("--valuecolumnidx") ) {
        valueColumnIdx = Integer.parseInt(argVal);
      } else if ( argName.equals("--cutoff") ) {
        cutoff = Double.parseDouble(argVal);
      } else if ( argName.equals("--taxid") ) {
        taxId = Integer.parseInt(argVal);        
      } else if ( argName.equals("--whatisbetter") ) {
        whatIsBetter = argVal;
      } else if ( argName.equals("--measure") ) {
        measure = argVal.toUpperCase();
      } else if ( argName.equals("--stdout") ) {
        stdout = true;
      }
    }
    
    GeneOntology myGo = new GeneOntology("go_201410-termdb.obo-xml.gz");
    
    GoTermSizes goSizes = null;
    try {
      goSizes = GoTermSizes.deserializeFromResources(goaPackageName, taxId);
    } catch (Exception e) {
      System.out.println("You have entered an unsupported TaxID. If Revigo doesn't directly support your organism, it's usually safe to use '0',"
              + " meaning the entire Uniprot will be used to determine GO Term frequencies. Aborting.");
      return;
    }
    
    goIdColumnIdx--;  // convert to 0-based column indices
    valueColumnIdx--;
 
    File inDirOrFile = new File(inputDir);
    File[] loopThrough;
    if ( inDirOrFile.isDirectory() ) {
      loopThrough = inDirOrFile.listFiles();
    } else {
      loopThrough = new File[1];
      loopThrough[0] = inDirOrFile; 
    }
    
    for ( File aFile : loopThrough ) {
      
      if ( aFile.getName().endsWith("_revigoOut.txt") ) {
        continue; // don't process Revigo output files
      }
      BufferedReader bufRdr = new BufferedReader(new FileReader(aFile));
      String line; 
      // load whole file into the userInput String
      List<String> userInputLines = new ArrayList<String>();
      while ((line = bufRdr.readLine()) != null) {
        line=line.trim();
        if ( line.isEmpty() || line.startsWith("%") || line.startsWith("#") ) {
          continue;
        }
        userInputLines.add(line);
      }
      
      GoTermList mfTerms = new GoTermList(), ccTerms = new GoTermList(),
              bpTerms = new GoTermList(), allTerms = new GoTermList();
      
      GoTermProperties termProps = new GoTermProperties(myGo);
      int termsWithValsCount = 0;

      List<String> warnings = new ArrayList<String>();
      int minNumColsPerGoTerm = Integer.MAX_VALUE;
      
      for (String curLine : userInputLines) {

        String[] cols = curLine.split("\t");
        int curGoId;
        double curVal;
        Integer curSize;
        //out.println(Arrays.toString(cols)+"<br>");

        cols[goIdColumnIdx] = cols[goIdColumnIdx].trim().replaceFirst("GO:", "").replaceFirst("GO", "");
        try {
          curGoId = Integer.parseInt(cols[goIdColumnIdx]);
        } catch (NumberFormatException e) {
          warnings.add("Could not parse GO ID from line: '" + curLine +
                  "'. Line will be skipped.");
          continue;
        }

        GOTerm curGoTerm = myGo.get(curGoId);
        if (curGoTerm == null) {
          warnings.add("Go term " + curGoId + " was not found in the "
                  + "current version of the GeneOntology, dated "
                  + myGo.getDate() + ". GO term will be skipped.");
          continue;
        }

        if (cols.length == 1) {
          curVal = Double.NaN;
        } else {
          try {
            curVal = Double.parseDouble(cols[valueColumnIdx].trim());
            termsWithValsCount++;
          } catch (NumberFormatException e) {
            warnings.add("Failed to parse number following GO ID in line: '"
                    + curLine + "'. This GO term will be treated as not having any value associated to it.");
            curVal = Double.NaN;
          }
        }
        
        // number of columns with user-supplied numeric data, in addition to the GO term
        // minNumColsPerGoTerm = Math.min(minNumColsPerGoTerm, cols.length-1); 
        // currently not supported in Batch Revigo, and always set to 1
        minNumColsPerGoTerm = 1;

        if (curGoTerm.isMolecularFunction()) {
          mfTerms.add(curGoTerm);
        } else if (curGoTerm.isCellularComponent()) {
          ccTerms.add(curGoTerm);
        } else if (curGoTerm.isBiologicalProcess()) {
          bpTerms.add(curGoTerm);
        }
        allTerms.add(curGoTerm);

        // out.println(curGoTerm.getId() + "<br/>");

        // add attributes to terms! (1) size and (2) value

        // if size not defined, attempts to guess using heuristics
        curSize = goSizes.getSizeGuessIfUnknown(curGoId, myGo); 

        termProps.put(curGoTerm, TermListVisualizer.log10SizeString, Math.log10(curSize));
        double transformedVal = Double.NaN;
        if ( whatIsBetter.equals("p_value") ) {
          if (curVal < 0) {
            throw new IllegalArgumentException("You can't specify a negative " +
                    "p-value. Go back to the input page and correct the error.");
          }
          curVal = Math.max(1e-300, curVal);
          transformedVal = -Math.log10(curVal);
          curVal = Math.log10(curVal);
        } else {
          if ( whatIsBetter.equals("higher") ) {
            transformedVal = curVal;
          } else if ( whatIsBetter.equals("lower") ) {
            transformedVal = -curVal;
          } else if ( whatIsBetter.equals("absolute") ) {
            transformedVal = Math.abs(curVal);
          } else if ( whatIsBetter.equals("abs_log") ) {
            if (curVal <= 0) {
              throw new IllegalArgumentException("Please provide only" +
                      " positive values if you" +
                      " choose the log-transform. Go back to the " +
                      "input page to correct the error.");
            }
            transformedVal = Math.abs( Math.log(curVal)/Math.log(2) );
            curVal = Math.log(curVal)/Math.log(2);
          }
        }

        if ( curVal != Double.NaN ) {
          termProps.put(curGoTerm, "value", curVal);
          termProps.put(curGoTerm, "transformedValue", transformedVal);
        }
        termProps.put(curGoTerm, "pinned", 0.0);

      } // line by line in user input file  
       
      // DEBUG -- set manually to "no additional data" - remove afterwards if needed
      // minNumColsPerGoTerm = 1;
      
      // check if user has consistently provided more than 1 numeric value per GO term
      // if yes, save it for later, and give it to table/graph which will be output later
      Map<GOTerm, double[]> allValsForTerms = null;
      // BatchREVIGO currently does not support this, minNumColsPerGoTerm always == 1
      if (minNumColsPerGoTerm > 1) {  
        allValsForTerms = new LinkedHashMap<GOTerm, double[]>();
        for (String curLine : userInputLines) {

          curLine = curLine.trim();
          //out.println(curLine+"<br>");
          if (curLine.length() == 0 || curLine.startsWith("%")) {
            continue;
          }
          String[] cols = curLine.split("\t");

          int curGoId;
          double curVal;
          cols[goIdColumnIdx] = cols[goIdColumnIdx].replaceFirst("GO:", "").replaceFirst("GO", "");
          try {
            curGoId = Integer.parseInt(cols[goIdColumnIdx]);
          } catch (NumberFormatException e) {
            continue;
          }
          GOTerm curGoTerm = myGo.get(curGoId);
          if ( ! allTerms.contains(curGoTerm) ) {
            continue;
          }
          double[] curArr = new double[cols.length-1];
          for (int c=0; c<cols.length; c++) {  // all columns except the one with the GO term
            if ( c==goIdColumnIdx )
              continue;
            try {
              curArr[c - 1] = Double.parseDouble(cols[c]);
            } catch (NumberFormatException e) {
              curArr[c - 1] = Double.NaN;
            }
          }
          allValsForTerms.put(curGoTerm, curArr);
      
        } // -------- line by line in input data --------
        
      } // ---- check if additional values were provided consistently
      
    
      boolean allTermsWithoutValue = false;
      
      if ( termsWithValsCount == 0 ) {
        allTermsWithoutValue = true;
      } else if ( termsWithValsCount < allTerms.size() ) {
        warnings.add("You have provided a numeric value only for some of the "
                + "GO terms ("+ termsWithValsCount +" out of "+allTerms.size()+"); "
                + "while REVIGO can still function, it is possible that you have"
                + " an error in your input data.");
      }

      String hugeListWarningTemplate = "You have provided a huge list of GO terms, containing "
                + "%d terms in the %s namespace, and only approx. "
                + TermListVisualizer.maxAllowedGoListSize +" of these terms will"
                + " be shown in the table below, or be available in the exported text table."
                + " REVIGO keeps the terms with the strongest p-values (or enrichements, depending what"
                + " you had specified on the input form). Consider filtering the"
                + " list by an external criterion (e.g. strength of enrichment) before"
                + " submitting to REVIGO.";
      if ( bpTerms.size() > TermListVisualizer.maxAllowedGoListSize  ) {
        warnings.add( String.format(hugeListWarningTemplate, bpTerms.size(), "BIOLOGICAL_PROCESS") );
      }
      if ( mfTerms.size() > TermListVisualizer.maxAllowedGoListSize  ) {
        warnings.add( String.format(hugeListWarningTemplate, mfTerms.size(), "MOLECULAR_FUNCTION") );
      }
      if ( ccTerms.size() > TermListVisualizer.maxAllowedGoListSize  ) {
        warnings.add( String.format(hugeListWarningTemplate, ccTerms.size(), "CELLULAR_COMPONENT") );
      }
      
      if ( warnings != null && warnings.size() > 0 ) {
        System.err.println("Warnings for file " + aFile );
        for ( String w : warnings ) {
          System.err.println(w);
        }
      }
      
      // make Visualizers for each namespace separately (and one for all together)
      Map<GoNamespace, TermListVisualizer> myVis = new HashMap<GoNamespace, TermListVisualizer>();
      myVis.put(GoNamespace.BIOLOGICAL_PROCESS,
              new TermListVisualizer(bpTerms, myGo, goSizes));
      myVis.put(GoNamespace.MOLECULAR_FUNCTION,
              new TermListVisualizer(mfTerms, myGo, goSizes));
      myVis.put(GoNamespace.CELLULAR_COMPONENT,
              new TermListVisualizer(ccTerms, myGo, goSizes));
      myVis.put(GoNamespace.MIXED_NAMESPACE,
              new TermListVisualizer(allTerms, myGo, goSizes));

      Integer minNumValsPerTerm = minNumColsPerGoTerm;
      Map<GOTerm, double[]> valsForTerms = allValsForTerms;

      PrintWriter out;
      if ( !stdout ) {
        // where to send the output.. to the same directory, to a file with the same name
        // with _revigoed.txt added to the end
        out = new PrintWriter( new FileWriter( aFile.getAbsolutePath()+"_revigoOut.txt" ) );
      } else {
        out = new PrintWriter( System.out );
      }
      
      if ( stdout ) {  // prints a spacer and an extra header line between different tables
        out.print("\n");
        out.print("%processing file: " + aFile.getName() + "\n");
      }
      
      // header
      out.print("%namespace\tterm_ID\tdescription\tfrequency\tplot_X\tplot_Y\tplot_size\t");
      if ( ! allTermsWithoutValue ) {
        out.print( ( whatIsBetter.equals("p_value") ? "log10_p-value" : ( whatIsBetter.equals("abs_log") ? "log2_value" : "value" ) ) + "\t" );
      }
      if (minNumValsPerTerm > 1) {  // this is currenty fixed at 1 for BatchREVIGO
        for (int c=1; c<minNumValsPerTerm; c++) {
          out.print("userVal_" + (c+1) + "\t" );
        }
      }
      out.println("uniqueness\tdispensability\trepresentative\teliminated");
      
      SemanticSimilarityScore scoreToUse;
      try {
        scoreToUse = SemanticSimilarityScore.valueOf(measure);
      } catch (Exception e) {
        System.out.println("You have entered an unsupported semantic similarity score. Aborting.");
        return;
      }
      
      for ( GoNamespace curNamespace : GoNamespace.values() ) {
        
        if (curNamespace == GoNamespace.MIXED_NAMESPACE || myVis.get(curNamespace).termList.size() == 0) {
          continue;
        }

        /* Used to store x,y,dispensability & co into properties. */
        myVis.get(curNamespace).makeOntologram(termProps, "transformedValue",
            cutoff,
            scoreToUse, goSizes );

        GoTermList terms = myVis.get(curNamespace).termList;
        terms.findClustersAndSortByThem(termProps, cutoff);

        for ( int i = 0; i<terms.size(); i++ ) { // row by row
          
          GOTerm aGoTerm = terms.get(i);
          boolean isTermEliminated = termProps.get(aGoTerm, "dispensability") > cutoff;
          out.print( aGoTerm.getNamespace().toString().replace(' ', '_') + "\t" );
          out.print( aGoTerm.getFormattedId() + "\t" );
          out.print( ""+aGoTerm.getName().replace(' ', '_')+"\t" );
          out.print( String.format( Locale.US, "%.3f%%\t", goSizes.getFreqGuessIfUnknown(aGoTerm.getId(),myGo)*100.0 ) );
          out.print( termProps.get(aGoTerm, "PC_1")==null ? "null\t" : String.format( Locale.US, "%6.3f\t", termProps.get(aGoTerm, "PC_1")) );
          out.print( termProps.get(aGoTerm, "PC_2")==null ? "null\t" : String.format( Locale.US, "%6.3f\t", termProps.get(aGoTerm, "PC_2")) );
          out.print( String.format( Locale.US, "%.3f\t", Math.log10( Math.max(1,goSizes.getSizeGuessIfUnknown(aGoTerm.getId(),myGo)) + 1 )) );
          if ( ! allTermsWithoutValue ) {
            out.print(String.format( Locale.US, "%.4f\t", termProps.get(aGoTerm, "value" ) ));
          }
          if (minNumValsPerTerm > 1) {
            double[] valsForCurTerm = valsForTerms.get(aGoTerm);
            if ( valsForCurTerm == null )
              throw new IllegalArgumentException("BatchREVIGO currently doesn't support >1 user-supplied numeric value (p-val, enrichment...) per GO term.");
            for (int c=1; c<minNumValsPerTerm; c++) {
              out.print( String.format( Locale.US, "%.3f\t", valsForCurTerm[c] ) );
            }
          }
          out.print( String.format( Locale.US, "%.3f\t", termProps.get(aGoTerm, "uniqueness") ) );
          out.print( String.format( Locale.US, "%.3f\t", termProps.get(aGoTerm, "dispensability") ) );
          out.printf( "GO:%07d\t", termProps.get(aGoTerm, "representative").intValue() );
          out.println( isTermEliminated ? "1" : "0");
          
        } // end row by row loop
        
      } // namespace by namespace loop     
      
      out.close();
      
    } // ==== END BIG LOOP ====  file by file
    
    
  }
  
  
  
}
